1 Install Packages

# # CRAN
# 
# install.packages("vegan")
# install.packages("Polychrome")
# install.packages("dendextend") 
# install.packages("ggplotify")
# install.packages("parallelMap")
# install.packages("caret")
# install.packages("randomForest")
# 
# # Bioconductor
# 
# install.packages("BiocManager")
# BiocManager::install("phyloseq")
# BiocManager::install("ANCOMBC")
# BiocManager::install("mixOmics")
# 
# # GitHub
# 
# install.packages("devtools")
# devtools::install_github("jbisanz/qiime2R")
library(dplyr)
library(ggplot2)

2 Load data to Phyloseq

Data is ASVs

# Required library
library(qiime2R) 
library(phyloseq) 
library(vegan)

# Source files:
feature_table_qza <- "feature_table.qza"
rooted_tree_qza <- "rooted_tree.qza"
taxonomy_qza <- "taxonomy.qza"
metadata_tsv <- "samples.txt"

# Read data
data.phy <- qza_to_phyloseq(
    features = feature_table_qza,
    tree = rooted_tree_qza,
    taxonomy = taxonomy_qza,
    metadata = metadata_tsv
    )
# Cleanup
# Removal of not needed objects, packages and cleaning the RAM
rm(feature_table_qza, rooted_tree_qza, taxonomy_qza, metadata_tsv) # remove unnecessary objects
detach("package:qiime2R", unload=TRUE) # unload qiime2R package
gc() # free unused R memory
##            used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells  5098515 272.3    8456080 451.7         NA  8456080 451.7
## Vcells 29988918 228.8   88611965 676.1      16384 91949714 701.6

3 Filtering and merging

3.1 Filter Industrial and Agriculture samples

# Keep only samples with Industrial or Agricultural prior use
table(sample_data(data.phy)$Former_landuse)
## 
##      Agriculture Ancient_woodland       Industrial        Rewilding 
##              150               20              150               10
IndAgri <- subset_samples(data.phy, 
                          Former_landuse %in% c("Industrial","Agriculture"))
IndAgri 
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 60355 taxa and 300 samples ]
## sample_data() Sample Data:       [ 300 samples by 11 sample variables ]
## tax_table()   Taxonomy Table:    [ 60355 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 60355 tips and 60236 internal nodes ]

3.2 Filter bacterial Kingdom

table(as.data.frame(tax_table(data.phy))$Kingdom)
## 
##  d__Archaea d__Bacteria  Unassigned 
##         161       60173          21
# Keep only Bacteria 
BacIndAgri <- subset_taxa(IndAgri, Kingdom == "d__Bacteria")
BacIndAgri # The OTU count dropped from 60,355 to 60,173
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 60173 taxa and 300 samples ]
## sample_data() Sample Data:       [ 300 samples by 11 sample variables ]
## tax_table()   Taxonomy Table:    [ 60173 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 60173 tips and 60054 internal nodes ]

3.3 Merging and Rarefaction

plot_rarefaction_curve <- function(phy.obj, taxa_level) {
  # sample depth Curve before Rarefaction
  # Extract the ASV table from the phyloseq object
  otu_table <- otu_table(phy.obj)
  # Transpose the table 
  if (taxa_are_rows(otu_table)) {
    otu_table <- t(otu_table)
  }
  # Convert to a matrix
  otu_matrix <- as(otu_table, "matrix")

  # Generate rarefaction curve to decide sample depth 
  rare_curve <- rarecurve(otu_matrix, step = 100, cex = 0.5, col = "blue", label = FALSE,
            xlab = "Sample Size", ylab = taxa_level)
}

3.4 1. Family level merging

# Agglomerate Bacteria to family level
BacIndAgriFamily <- tax_glom(BacIndAgri, taxrank="Family")
BacIndAgriFamily
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 769 taxa and 300 samples ]
## sample_data() Sample Data:       [ 300 samples by 11 sample variables ]
## tax_table()   Taxonomy Table:    [ 769 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 769 tips and 768 internal nodes ]

3.4.1 Normalization (Rarefaction)

# Check the sample depth for deciding the rarefaction threshold
plot_rarefaction_curve(BacIndAgriFamily, "ASVs (Family)")

3.5 2. Genera level merging

# Agglomerate Bacteria to genus level
BacIndAgriGenera <- tax_glom(BacIndAgri, taxrank="Genus")
BacIndAgriGenera
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1465 taxa and 300 samples ]
## sample_data() Sample Data:       [ 300 samples by 11 sample variables ]
## tax_table()   Taxonomy Table:    [ 1465 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1465 tips and 1464 internal nodes ]

3.5.1 Normalization (Rarefaction)

# Check the sample depth for deciding the rarefaction threshold
plot_rarefaction_curve(BacIndAgriGenera, "ASVs (Genus)")

3.6 3. Species level merging

3.6.1 Normalization (Rarefaction)

# Check the sample depth for deciding the rarefaction threshold
plot_rarefaction_curve(BacIndAgri, "ASVs (Species)")

3.7 Applying Rarefaction on Genera merged data

# Rarefy the phyloseq object to an even depth of 10000 sequences per sample
BacIndAgriGenera.rarefied <- rarefy_even_depth(BacIndAgriGenera, 
                                                 sample.size = 10000, 
                                                 rngseed = 123, replace = FALSE)
  • Rarefation done on Genera level, for multiple reasons: 1. not all ASVs are classified up to species 2. The 16srRNA sequencing achieves better genus-level resolution

  • The decision of sample depth threshold (10000) to use for normalization was mainly based on visual assessment of two things: sample reaching the plateau of ASVs count and avoiding losing any samples.

3.8 Filtering low quality ASVs (at Genera level)

# (Function) to create and save a histogram
save_histogram <- function(data, xlab, ylab, main, filename) {
  # Display the histogram
  hist(data, xlab = xlab, ylab = ylab, main = main)
  
  # Save the histogram to a file
  dev.copy(png, filename)
  dev.off()
}
# Convert the OTU table to matrix
BacIndAgriGenera.numeric <- as.numeric(otu_table(BacIndAgriGenera))
BacIndAgriGenera.mtx <- matrix(BacIndAgriGenera.numeric, nrow=nrow(otu_table(BacIndAgriGenera)))

# Getting the ASV that has more than 0 count in each sample THRESHOLD
BacIndAgriGenera.count_filter <- BacIndAgriGenera.mtx > 0

# Sum the number of ASV for each sample that has > count threshold 
# The result will show the number of ASV (frequency, y-axis) > count threshold in each sample (x-axis)
BacIndAgriGenera.count_filter.sum <- apply(BacIndAgriGenera.count_filter, 1, sum)

# Display and save the first histogram
save_histogram(BacIndAgriGenera.count_filter.sum, xlab = "Samples", ylab = "ASVs count (Genera)",
               main = "Frequency of ASVs with count > 0 across samples", 
               filename = "results/filter_low_quality_ASVs/ASVsGenera_count_histogram.png")

## quartz_off_screen 
##                 2
# Plot histogram for sample interval (0-10)
BacIndAgriGenera.count_filter.sum.1_10.interval <- BacIndAgriGenera.count_filter.sum[BacIndAgriGenera.count_filter.sum < 10]

# Display and save the second histogram
save_histogram(BacIndAgriGenera.count_filter.sum.1_10.interval, xlab = "Samples", ylab = "ASVs count (Genera)",
               main = "Frequency of ASVs with count > 0 in 10 samples or less", 
               filename = "results/filter_low_quality_ASVs/ASVsGenera_count_histogram_1_10_interval.png")

## quartz_off_screen 
##                 2

Interpretation: Let’s say for the second plot we have ~ 150 ASVs that are present in two samples.

all_taxa.num <- length(BacIndAgriGenera.count_filter.sum)
taxa_less10samples <- length(BacIndAgriGenera.count_filter.sum.1_10.interval)
lost_taxa_perc <- round((taxa_less10samples / all_taxa.num) * 100, 2)

print(paste("Out of", all_taxa.num, "the number of lost ASVs (Genera) after filtering taxa with abundance > 0 counts in less than 10 samples", taxa_less10samples, "Which is", lost_taxa_perc, "%"))
## [1] "Out of 1465 the number of lost ASVs (Genera) after filtering taxa with abundance > 0 counts in less than 10 samples 766 Which is 52.29 %"

From count filtering histogram and stats, half of ASVs (Genera) will be lost. Therefore, for further analysis the data will be kept without filtering unless it’s required for a statistical method to be used.

4 Visualizing Microbial diversity

library(ggplot2)
library(plotly)
library(RColorBrewer)

4.1 Bar plot

4.1.1 Former land use

# Aggregate data at Phylum level
BacIndAgriPhylum.rarefied <- tax_glom(BacIndAgriGenera.rarefied, taxrank = "Phylum")

4.2 Check data types

# Make sure pH / Conductivity / Age are numerical 
sample_data(BacIndAgriPhylum.rarefied)$pH <- as.numeric(as.character(sample_data(BacIndAgriPhylum.rarefied)$pH))
sample_data(BacIndAgriPhylum.rarefied)$Conductivity <- as.numeric(as.character(sample_data(BacIndAgriPhylum.rarefied)$Conductivity))
sample_data(BacIndAgriPhylum.rarefied)$Woodland_age <- as.numeric(as.character(sample_data(BacIndAgriPhylum.rarefied)$Woodland_age))
env_factor_bar_plot <- function(phyloseq.obj, env_factor, save_path, levels=c()) { 
  # Merge samples by environmental factor (categorical one)
  phyloseq.obj.merged <- merge_samples(phyloseq.obj, env_factor)
  
  # Transform counts to relative abundances
  phyloseq.obj.merged <- transform_sample_counts(phyloseq.obj.merged, function(x) x / sum(x))
  
  # Extract the data for plotting
  phyloseq.obj.merged.df <- psmelt(phyloseq.obj.merged)
  
  
  # Get counts per phylum
  Phyla.df <- phyloseq.obj.merged.df %>% 
    group_by(Phylum) %>% 
    summarise(Count = sum(Abundance))
  
  # Select the cut-off for the Phylum taxa (e.g. 1% of total count)
  cutoff <- 0.01 * sum(phyloseq.obj.merged.df$Abundance)
  
  # Select low-abundant Phyla (with total counts below the cutoff)
  lowAbundant <- Phyla.df[Phyla.df$Count <= cutoff,]$Phylum
  
  # Substitute Phylum names to "<1%" for the low-abundant phyla
  phyloseq.obj.merged.df[phyloseq.obj.merged.df$Phylum %in% lowAbundant,]$Phylum <- '<1%'
  
  # Ensure env_factor levels are ordered correctly in the melted data frame
  if (length(levels) > 0){
    phyloseq.obj.merged.df$Sample <- factor(phyloseq.obj.merged.df$Sample, levels = levels)
  }
  
  # Define high contrast palette for bar plot
  n_colors <- length(unique(phyloseq.obj.merged.df$Phylum))
  high_contrast_palette <- c(brewer.pal(min(9, n_colors), "Set1"), brewer.pal(max(0, n_colors - 9), "Dark2"))
  
  # Create bar plot 
  barplot <- ggplot(phyloseq.obj.merged.df, aes(x = Sample, y = Abundance, fill = Phylum, text = paste("Phylum:", Phylum, "<br>", env_factor, ":", Sample, "<br>Relative Abundance:", Abundance))) +
    geom_bar(stat = "identity", position = "stack") +
    scale_fill_manual(values = high_contrast_palette) +
    theme(axis.text.x = element_text(angle = 0, hjust = 1)) +
    labs(title = paste("Taxa Abundance by", env_factor), x = env_factor, y = "Relative Abundance")
  
  # Convert to plotly 
  barplot.plotly <- ggplotly(barplot, tooltip = "text")
  barplot.plotly
  
  # Save the plot as an HTML file
  htmlwidgets::saveWidget(barplot.plotly, save_path)
  return(barplot.plotly)
} 
barplot.FLU <- env_factor_bar_plot(phyloseq.obj = BacIndAgriPhylum.rarefied, env_factor = "Former_landuse",
                    save_path = "results/plotly/Taxonomy_Phyla_landuse_normalized.html")
barplot.FLU

4.2.1 pH

pHValues <- as.numeric(sample_data(BacIndAgriGenera)$pH)
## Warning: NAs introduced by coercion
pHValues.removed_nan <- na.omit(pHValues)
print("")
## [1] ""
print("max pH")
## [1] "max pH"
max(pHValues.removed_nan)
## [1] 7.98
print("Min pH")
## [1] "Min pH"
min(pHValues.removed_nan)
## [1] 3.78
print("Median pH")
## [1] "Median pH"
median(pHValues.removed_nan)
## [1] 5.54
# Make sure pH column is numeric
sample_data(BacIndAgriPhylum.rarefied)$pH <- as.numeric(as.character(sample_data(BacIndAgriPhylum.rarefied)$pH))

# Categorizing pH 
sample_data(BacIndAgriPhylum.rarefied)$pH_category <- factor(
  cut(sample_data(BacIndAgriPhylum.rarefied)$pH,
      breaks = c(3.78, 5, 6, 7.98),
      labels = c("<5", "5-6", ">6")),
  levels = c("<5", "5-6", ">6")
)
pH_less_5 <- as.numeric(sample_data(BacIndAgriPhylum.rarefied)$pH_category == "<5")
pH_less_5 <- na.omit(pH_less_5)
print(paste("Number of samples with to pH < 5 =", sum(pH_less_5)))
## [1] "Number of samples with to pH < 5 = 50"
pH_5_6 <- as.numeric(sample_data(BacIndAgriPhylum.rarefied)$pH_category == "5-6")
pH_5_6 <- na.omit(pH_5_6)
print(paste("Number of samples with to pH 5-6 =", sum(pH_5_6)))
## [1] "Number of samples with to pH 5-6 = 145"
pH_more_6 <- as.numeric(sample_data(BacIndAgriPhylum.rarefied)$pH_category == ">6")
pH_more_6 <- na.omit(pH_more_6)
print(paste("Number of samples with pH >6 =", sum(pH_more_6)))
## [1] "Number of samples with pH >6 = 95"
barplot.pH <- env_factor_bar_plot(phyloseq.obj = BacIndAgriPhylum.rarefied, env_factor = "pH_category",
                    save_path = "results/plotly/Taxonomy_Phyla_pH_normalized.html",
                    levels = c("<5", "5-6", ">6"))
barplot.pH

4.2.2 Conductivity

conductivityValues <- as.numeric(sample_data(BacIndAgriGenera)$Conductivity)
## Warning: NAs introduced by coercion
# Calculate the quartiles
quartiles <- quantile(conductivityValues, probs = c(0.25, 0.5, 0.75), na.rm = TRUE)
quartiles
## 25% 50% 75% 
##  39  64 134
# Make sure conductivity column is numeric
sample_data(BacIndAgriPhylum.rarefied)$Conductivity <- as.numeric(as.character(sample_data(BacIndAgriPhylum.rarefied)$Conductivity))

# Calculate the quartiles for conductivity
quartiles <- quantile(sample_data(BacIndAgriPhylum.rarefied)$Conductivity, probs = c(0.25, 0.5, 0.75), na.rm = TRUE)

# Categories based on quartiles
sample_data(BacIndAgriPhylum.rarefied)$Conductivity_category <-
  cut(sample_data(BacIndAgriPhylum.rarefied)$Conductivity,
      breaks = c(-Inf, quartiles, Inf),
      labels = c("Q1: <0.25", "Q2: 0.25-0.5", "Q3: 0.5-0.75", "Q4: 0.75<"),
      include.lowest = TRUE)
barplot.conductivity <- env_factor_bar_plot(phyloseq.obj = BacIndAgriPhylum.rarefied, env_factor = "Conductivity_category",
                    save_path = "results/plotly/Taxonomy_Phyla_Conductivity_normalized.html")
barplot.conductivity

4.3 Line Plot

4.3.1 FLU vs pH

# Make sure pH column is numeric
sample_data(BacIndAgri)$pH <- as.numeric(as.character(sample_data(BacIndAgri)$pH))

# Make Former_landuse column as factor
sample_data(BacIndAgri)$Former_landuse <- as.factor(sample_data(BacIndAgri)$Former_landuse)

# Extract sample data into a df
FLU_pH.df <- data.frame(sample_data(BacIndAgri))

FLU_pH.plot <- plot_ly()

# Compute density estimates for each Former_landuse level
for(level in levels(FLU_pH.df$Former_landuse)) {
  data_subset <- FLU_pH.df[FLU_pH.df$Former_landuse == level, ]
  
  dens <- density(data_subset$pH, na.rm = TRUE)  # Compute density
  # Add density trace
  FLU_pH.plot <- FLU_pH.plot %>% add_trace(
    x = dens$x,
    y = dens$y,
    type = 'scatter',
    mode = 'lines',
    name = level,
    opacity = 0.6
  )
}

# Customize the layout
FLU_pH.plot <- FLU_pH.plot %>% layout(
  title = "Kernel Density Estimate of pH by Former Landuse",
  xaxis = list(title = "pH"),
  yaxis = list(title = "Density")
)

FLU_pH.plot
# Save the plot as an HTML file
htmlwidgets::saveWidget(FLU_pH.plot, "results/plotly/pH_vs_FLU.html")

4.4 Correlation

Scatter Plot ### pH vs Conductivity

library(broom)
# Make sure conductivity column is numeric
sample_data(BacIndAgri)$Conductivity <- as.numeric(as.character(sample_data(BacIndAgri)$Conductivity))
## Warning: NAs introduced by coercion
# Extract sample data 
Cond_pH.df <- data.frame(sample_data(BacIndAgri))

# linear regression 
model <- lm(Conductivity ~ pH, data = Cond_pH.df)

# Extract the coefficients
summary_lm <- summary(model) # Benjamini-Hochberg adjusted p-value by default
print(summary_lm)
## 
## Call:
## lm(formula = Conductivity ~ pH, data = Cond_pH.df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -85.67 -37.37 -11.57  29.63 237.12 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -96.885     20.258  -4.783 2.74e-06 ***
## pH            31.986      3.474   9.207  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 56.62 on 293 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.2244, Adjusted R-squared:  0.2218 
## F-statistic: 84.77 on 1 and 293 DF,  p-value: < 2.2e-16
# Extract R-squared and slope
r_squared <- summary_lm$r.squared
slope <- summary_lm$coefficients[2, 1]

# Create the regression plot
Cond_pH.plot <- ggplot(Cond_pH.df, aes(x = pH, y = Conductivity, color = Former_landuse)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(title = "Regression of Conductivity on pH by Former Land Use",
       subtitle = paste("R-squared =", round(r_squared, 2), 
                          "Slope =", round(slope, 2)),
       x = "pH",
       y = "Conductivity") +
  theme_minimal()

print(Cond_pH.plot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Convert to plotly 
Cond_pH.plot_plotly <- ggplotly(Cond_pH.plot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_smooth()`).
print(Cond_pH.plot_plotly)

# Save the plot as an HTML file
htmlwidgets::saveWidget(Cond_pH.plot_plotly, "results/plotly/Regression_pH_Conductivity_by_LandUse.html")

Adjusted R-squared 0.2218 (very low) shows that the model failed to explain the real data points and the results are not reliable. For the results it tells that for one unit increase of pH the conductivity increases by 31.986, and theoretically at pH 0 the conductivity is -96.885

4.4.1 Woodland age vs taxa abundance

library(dplyr)
# Function to compute the correlation between microbial abundance and environmental factor
abundance_env_factor_corr <- function(phyloseq.obj, env_factor){
  # Create abundance df
  phyloseq.obj.abund <- psmelt(phyloseq.obj)
  
  # Ensure environmental factors are numerical
  phyloseq.obj.abund$Woodland_age <- as.numeric(as.character(phyloseq.obj.abund$Woodland_age))
  phyloseq.obj.abund$Conductivity <- as.numeric(as.character(phyloseq.obj.abund$Conductivity))
  phyloseq.obj.abund$pH <- as.numeric(as.character(phyloseq.obj.abund$pH))
  
  # Aggregate abundance by Site_code (as samples of each site share the same value of pH, Conductivity, age) and Phylum, summing ASVs
  phyloseq.obj.abund.summary <- phyloseq.obj.abund %>%
    group_by(Site_code, Phylum) %>%
    summarise(Abundance = mean(Abundance), .groups = 'drop')
  
  # Align Sample IDs between phyloseq sample data and the summary data
  sample_data_df <- as.data.frame(sample_data(phyloseq.obj))
  sample_data_df$Sample <- rownames(sample_data_df)
  
  # Join the summary data with the sample metadata
  phyloseq.obj.abund.summary <- left_join(phyloseq.obj.abund.summary, sample_data_df, by = "Site_code")
  
  # Correlation function
  compute_correlation <- function(df) {
    cor_result <- cor.test(df[[env_factor]], df$Abundance, method = "spearman")
    data.frame(
      corr = cor_result$estimate,
      p_value = cor_result$p.value
    )
  }
  
  # Compute correlation and p-value for each Phylum
  correlations <- phyloseq.obj.abund.summary %>%
    group_by(Phylum) %>%
    summarise(correlation = list(compute_correlation(cur_data())), .groups = 'drop')
  
  # Flatten the list column to extract correlation and p-value
  correlations <- correlations %>%
    mutate(
      correlation = sapply(correlations$correlation, function(x) x$corr),
      p_value = sapply(correlations$correlation, function(x) x$p_value)
    )
  # Calculate the adjusted p-value
  correlations$pAdjust <- p.adjust(correlations$p_value)
  
  # Filter taxa based on adjusted p-value <= 0.05
  significant_taxa <- correlations %>%
    filter(pAdjust <= 0.05)
  
  # Get top positively and negatively correlated taxa
  top_positive <- significant_taxa %>%
    arrange(desc(correlation)) %>%
    slice_head(n = 4)
  
  top_negative <- significant_taxa %>%
    arrange(correlation) %>%
    slice_head(n = 4)
  
  print("Top positively correlated taxa with Woodland age:")
  print(top_positive)
  
  print("Top negatively correlated taxa with Woodland age:")
  print(top_negative)
  
  return(list("abund.summary"=phyloseq.obj.abund.summary, 
          "top_positive"=top_positive, "top_negative"=top_negative))

}
age_abund_corr <- abundance_env_factor_corr(phyloseq.obj = BacIndAgriPhylum.rarefied, 
                                          env_factor = "Woodland_age")
## Warning in left_join(phyloseq.obj.abund.summary, sample_data_df, by = "Site_code"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 146 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
## Warning: There were 58 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `correlation = list(compute_correlation(cur_data()))`.
## ℹ In group 1: `Phylum = "Acidobacteriota"`.
## Caused by warning:
## ! `cur_data()` was deprecated in dplyr 1.1.0.
## ℹ Please use `pick()` instead.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 57 remaining warnings.
## [1] "Top positively correlated taxa with Woodland age:"
## # A tibble: 4 × 4
##   Phylum             correlation  p_value       pAdjust
##   <chr>                    <dbl>    <dbl>         <dbl>
## 1 CSP1-3                   0.370 3.70e-11 0.00000000211
## 2 Sumerlaeota              0.304 7.64e- 8 0.00000428   
## 3 Desulfobacterota_C       0.299 1.29e- 7 0.00000707   
## 4 KSB1                     0.283 6.44e- 7 0.0000348    
## [1] "Top negatively correlated taxa with Woodland age:"
## # A tibble: 4 × 4
##   Phylum              correlation    p_value   pAdjust
##   <chr>                     <dbl>      <dbl>     <dbl>
## 1 Firmicutes_E             -0.278 0.00000104 0.0000552
## 2 Firmicutes_A             -0.261 0.00000468 0.000234 
## 3 Firmicutes_B_370531      -0.222 0.000106   0.00454  
## 4 Planctomycetota          -0.199 0.000520   0.0218
top_pos_neg_plot <- function(top_positive, top_negative, abund.summary, env_factor) {
  # Select the top positively and negatively correlated taxa
  top_taxa <- rbind(top_positive[1, ], top_negative[1, ])
  
  # Extract data for these taxa
  top_taxa.abund.summary <- abund.summary %>%
    filter(Phylum %in% top_taxa$Phylum)
  
  
  # Scatter plots
  for (phylum in unique(top_taxa$Phylum)) {
    # Filter data for the current Phylum
    phylum_data <- top_taxa.abund.summary %>% filter(Phylum == phylum)
    
    # Fit a linear model
    lm_formula <- paste("Abundance", "~", env_factor)
    lm_formula <- as.formula(lm_formula)
    lm_model <- lm(lm_formula, data = phylum_data)
    summary_lm <- summary(lm_model)
    print("--------------------------------------")
    print(summary_lm)
    
    # Extract R-squared and slope
    r_squared <- summary_lm$r.squared
    slope <- summary_lm$coefficients[2, 1]
    
    # Create scatter plot
    p <- ggplot(phylum_data, aes_string(x = env_factor, y = "Abundance")) +
      geom_point(color = "blue") +
      geom_smooth(method = "lm", se = TRUE, color = "red") +
      labs(title = paste("Abundance of", phylum, "vs", env_factor),
           subtitle = paste("R-squared =", round(r_squared, 2), 
                            "Slope =", round(slope, 2)),
           x = env_factor,
           y = "Abundance") +
      theme_minimal()
    
    # Print the plot
    print(p)
  }
}
top_pos_neg_plot(top_positive = age_abund_corr$top_positive, top_negative = age_abund_corr$top_negative,
                 abund.summary = age_abund_corr$abund.summary, env_factor = "Woodland_age")
## [1] "--------------------------------------"
## 
## Call:
## lm(formula = lm_formula, data = phylum_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27042 -0.10884 -0.02256  0.03078  0.99233 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.150011   0.034941  -4.293 2.38e-05 ***
## Woodland_age  0.006275   0.000851   7.373 1.65e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2391 on 298 degrees of freedom
## Multiple R-squared:  0.1543, Adjusted R-squared:  0.1515 
## F-statistic: 54.37 on 1 and 298 DF,  p-value: 1.654e-12
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `geom_smooth()` using formula = 'y ~ x'
## [1] "--------------------------------------"
## 
## Call:
## lm(formula = lm_formula, data = phylum_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5679 -0.3770 -0.2071 -0.0433  4.7929 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.696519   0.117930   5.906 9.54e-09 ***
## Woodland_age -0.008039   0.002872  -2.799  0.00547 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.807 on 298 degrees of freedom
## Multiple R-squared:  0.02561,    Adjusted R-squared:  0.02234 
## F-statistic: 7.832 on 1 and 298 DF,  p-value: 0.005467
## `geom_smooth()` using formula = 'y ~ x'

4.4.2 pH vs taxa abundance

pH_abund_corr <- abundance_env_factor_corr(phyloseq.obj = BacIndAgriPhylum.rarefied, 
                                          env_factor = "pH")
## Warning in left_join(phyloseq.obj.abund.summary, sample_data_df, by = "Site_code"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 146 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
## Warning: There were 57 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `correlation = list(compute_correlation(cur_data()))`.
## ℹ In group 1: `Phylum = "Acidobacteriota"`.
## Caused by warning in `cor.test.default()`:
## ! Cannot compute exact p-value with ties
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 56 remaining warnings.
## [1] "Top positively correlated taxa with Woodland age:"
## # A tibble: 4 × 4
##   Phylum               correlation  p_value  pAdjust
##   <chr>                      <dbl>    <dbl>    <dbl>
## 1 Chloroflexota              0.699 1.34e-44 7.65e-43
## 2 Methylomirabilota          0.506 1.53e-20 8.10e-19
## 3 Myxococcota_A_437813       0.438 3.02e-15 1.51e-13
## 4 Myxococcota_A_473307       0.425 2.25e-14 1.10e-12
## [1] "Top negatively correlated taxa with Woodland age:"
## # A tibble: 4 × 4
##   Phylum          correlation  p_value  pAdjust
##   <chr>                 <dbl>    <dbl>    <dbl>
## 1 Eremiobacterota      -0.678 5.35e-41 2.99e-39
## 2 Planctomycetota      -0.675 1.23e-40 6.78e-39
## 3 Acidobacteriota      -0.559 1.17e-25 6.31e-24
## 4 Firmicutes_D         -0.505 1.58e-20 8.20e-19
top_pos_neg_plot(top_positive = pH_abund_corr$top_positive, top_negative = pH_abund_corr$top_negative,
                 abund.summary = pH_abund_corr$abund.summary, env_factor = "pH")
## [1] "--------------------------------------"
## 
## Call:
## lm(formula = lm_formula, data = phylum_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -273.56  -95.68   -8.36   69.29  497.47 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -200.772     54.059  -3.714 0.000244 ***
## pH           106.065      9.271  11.441  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 151.1 on 293 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.3088, Adjusted R-squared:  0.3064 
## F-statistic: 130.9 on 1 and 293 DF,  p-value: < 2.2e-16
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
## [1] "--------------------------------------"
## 
## Call:
## lm(formula = lm_formula, data = phylum_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.705  -8.128  -2.338   5.547  48.882 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  58.6198     4.1269   14.20   <2e-16 ***
## pH           -8.2371     0.7077  -11.64   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.53 on 293 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.3162, Adjusted R-squared:  0.3138 
## F-statistic: 135.5 on 1 and 293 DF,  p-value: < 2.2e-16
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 5 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).

4.4.3 Conductivity vs taxa abundance

conductivity_abund_corr <- abundance_env_factor_corr(phyloseq.obj = BacIndAgriPhylum.rarefied, 
                                                     env_factor = "Conductivity")
## Warning in left_join(phyloseq.obj.abund.summary, sample_data_df, by = "Site_code"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 146 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
## Warning: There were 57 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `correlation = list(compute_correlation(cur_data()))`.
## ℹ In group 1: `Phylum = "Acidobacteriota"`.
## Caused by warning in `cor.test.default()`:
## ! Cannot compute exact p-value with ties
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 56 remaining warnings.
## [1] "Top positively correlated taxa with Woodland age:"
## # A tibble: 4 × 4
##   Phylum              correlation  p_value  pAdjust
##   <chr>                     <dbl>    <dbl>    <dbl>
## 1 Actinobacteriota          0.699 1.69e-44 9.66e-43
## 2 Chloroflexota             0.546 2.66e-24 1.33e-22
## 3 Firmicutes_B_370539       0.454 2.16e-16 9.94e-15
## 4 Methylomirabilota         0.415 9.78e-14 4.40e-12
## [1] "Top negatively correlated taxa with Woodland age:"
## # A tibble: 4 × 4
##   Phylum          correlation  p_value  pAdjust
##   <chr>                 <dbl>    <dbl>    <dbl>
## 1 Acidobacteriota      -0.678 4.25e-41 2.38e-39
## 2 FCPU426              -0.677 5.91e-41 3.25e-39
## 3 Elusimicrobiota      -0.639 3.03e-35 1.64e-33
## 4 Armatimonadota       -0.632 2.47e-34 1.31e-32
top_pos_neg_plot(top_positive = conductivity_abund_corr$top_positive, 
                 top_negative = conductivity_abund_corr$top_negative,
                 abund.summary = conductivity_abund_corr$abund.summary, env_factor = "Conductivity")
## [1] "--------------------------------------"
## 
## Call:
## lm(formula = lm_formula, data = phylum_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2040.42  -602.79   -96.51   430.54  2617.57 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1432.4810    80.4577   17.80   <2e-16 ***
## Conductivity    9.3798     0.7438   12.61   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 818.5 on 293 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.3518, Adjusted R-squared:  0.3496 
## F-statistic:   159 on 1 and 293 DF,  p-value: < 2.2e-16
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
## [1] "--------------------------------------"
## 
## Call:
## lm(formula = lm_formula, data = phylum_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1269.72  -351.81   -18.32   318.38  1465.93 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2601.1759    52.4365   49.61   <2e-16 ***
## Conductivity   -5.5302     0.4848  -11.41   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 533.5 on 293 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.3075, Adjusted R-squared:  0.3052 
## F-statistic: 130.1 on 1 and 293 DF,  p-value: < 2.2e-16
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 5 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).

5 Alpha Diversity

5.1 plot

alpha_diversity.plot <- function(alpha_richness, phyloseq_obj, env_factor){
  # Apply Wilcoxon test 
  chao1_res <- wilcox.test(alpha_richness$Chao1~
                             sample_data(phyloseq_obj)[[env_factor]])
  chao1.p_value <- chao1_res$p.value
  shannon_res <- wilcox.test(alpha_richness$Shannon~
                               sample_data(phyloseq_obj)[[env_factor]])
  shannon.p_value <- shannon_res$p.value
  simpson_res <- wilcox.test(alpha_richness$Simpson~
                               sample_data(phyloseq_obj)[[env_factor]])
  simpson.p_value <- simpson_res$p.value
  
  plot_richness(phyloseq_obj, 
                measures=c("Chao1", "Shannon", "Simpson"),
                color=env_factor, x=env_factor) + 
    geom_boxplot() + ggtitle(paste("Alpha diversity vs", env_factor), 
                             subtitle=paste("Chao1 p-value:", chao1.p_value,
                                            "\nShannon p-value:", shannon.p_value,
                                            "\nSimpson p-value:", simpson.p_value))
}

5.1.1 by FLU

# Calculate Alpha diversity indices
alpha_richness = estimate_richness(
  BacIndAgriGenera.rarefied, measures = c("Chao1", "Shannon", "Simpson"))
alpha_diversity.plot(alpha_richness = alpha_richness, phyloseq_obj = BacIndAgriGenera.rarefied,
                     env_factor = "Former_landuse")

### by Region

alpha_diversity.plot(alpha_richness = alpha_richness, phyloseq_obj = BacIndAgriGenera.rarefied,
                     env_factor = "Region")

5.2 Variance

5.2.1 Shannon

shannon_values <- alpha_richness$Shannon
site_codes <- sample_data(BacIndAgriGenera.rarefied)$Site_code
FLU_values <- sample_data(BacIndAgriGenera.rarefied)$Former_landuse
shannon_site.df <- data.frame(Site_code = site_codes, Shannon = shannon_values, 
                              Former_landuse = FLU_values)
ggplot(shannon_site.df, aes(x = Site_code, y = Shannon, fill = Former_landuse)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Shannon Diversity by Site",
       x = "Site",
       y = "Shannon Index") +
  theme(axis.text.x = element_text(angle = 90, size=4)) 

# Calculate variance for sites replicas
variance_data.shannon <- shannon_site.df %>%
  group_by(Site_code) %>%
  summarise(CV = (sd(Shannon, na.rm = TRUE)/mean(Shannon, na.rm = TRUE))*100, Former_landuse = unique(Former_landuse))
ggplot(variance_data.shannon, aes(x = Site_code, y = CV, fill = Former_landuse)) +
  geom_bar(stat = "identity") + 
  theme_minimal() +
  labs(title = "Coefficient of Variation of Shannon Diversity by Site",
       x = "Site",
       y = "Coefficient of Variation") +
  theme(axis.text.x = element_text(angle = 90, size=4)) 

library(readxl)
sites_pH_EC <- read_excel("pH_EC_RestREco_8th_July_2024.xlsx")

5.2.2 pH

ggplot(sites_pH_EC, aes(x = Site_code, y = pH, fill = Former_landuse)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "pH Diversity by Site",
       x = "Site",
       y = "pH") +
  theme(axis.text.x = element_text(angle = 90, size = 6)) 

# Calculate variance for sites replicas
variance_data.pH <- sites_pH_EC %>%
  group_by(Site_code) %>%
  summarise(CV = (sd(pH, na.rm = TRUE)/mean(pH, na.rm = TRUE))*100, Former_landuse = unique(Former_landuse))
ggplot(variance_data.pH, aes(x = Site_code, y = CV, fill = Former_landuse)) +
  geom_bar(stat = "identity") + 
  theme_minimal() +
  labs(title = "Coefficient of Variation of pH by Site",
       x = "Site",
       y = "Coefficient of Variation") +
  theme(axis.text.x = element_text(angle = 90, size = 6)) 

5.2.3 Conductivity

ggplot(sites_pH_EC, aes(x = Site_code, y = Conductivity, fill = Former_landuse)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Conductivity Diversity by Site",
       x = "Site",
       y = "Conductivity") +
  theme(axis.text.x = element_text(angle = 90, size = 6)) 

# Calculate variance for sites replicas
variance_data.conductivity <- sites_pH_EC %>%
  group_by(Site_code) %>%
  summarise(CV = (sd(Conductivity, na.rm = TRUE)/mean(Conductivity, na.rm = TRUE))*100,
            Former_landuse = unique(Former_landuse))
ggplot(variance_data.conductivity, aes(x = Site_code, y = CV, fill = Former_landuse)) +
  geom_bar(stat = "identity") +  
  theme_minimal() +
  labs(title = "Coefficient of Variation of Conductivity by Site",
       x = "Site",
       y = "Coefficient of Variation") +
  theme(axis.text.x = element_text(angle = 90, size = 6)) 

5.3 Modelling (nested design)

#install.packages("remotes")
#remotes::install_github("mikemc/speedyseq")
library(speedyseq)
# Merge replicas for each site by median (non-parametric)
BacIndAgriGenera.rarefied.merged_replica <- merge_samples2(x = BacIndAgriGenera.rarefied, group = "Site_code", fun_otu = median)
# Sanity check of counts per sample (the sum should be same across all samples)
plot(colSums(otu_table(BacIndAgriGenera.rarefied.merged_replica)), 
     ylab="Count", xlab="Samples", xaxt="n", 
     main="Counts per sample after merging")

# Plot the ASVs count per site after median merging to decide a threshold for normalization
plot_rarefaction_curve(BacIndAgriGenera.rarefied.merged_replica, "ASVs (Genera)")

# Rarefy (normalize) to an even depth of 8000 sequences per sample
BacIndAgriGenera.rarefied2.merged_replica <- rarefy_even_depth(BacIndAgriGenera.rarefied.merged_replica, 
                                                 sample.size = 8000, 
                                                 rngseed = 123, replace = FALSE)
# Sanity check 
plot(colSums(otu_table(BacIndAgriGenera.rarefied2.merged_replica)), 
     ylab="Count", xlab="Samples", xaxt="n", 
     main="Counts per sample after merging and rarefaction")

# Calculate Shannon
richness.shannon = estimate_richness(
  BacIndAgriGenera.rarefied2.merged_replica, measures = c("Shannon"))
# Add shannon to the phyloseq obj
sample_data(BacIndAgriGenera.rarefied2.merged_replica)$Shannon <- richness.shannon$Shannon
# Get sample data into dataframe
BacIndAgriGenera.rarefied2.merged_replica.df <-
as(sample_data(BacIndAgriGenera.rarefied2.merged_replica), "data.frame")
reg_model.alphaDiv <- lm(Shannon~Region/Former_landuse, data=BacIndAgriGenera.rarefied2.merged_replica.df)
summary(reg_model.alphaDiv)
## 
## Call:
## lm(formula = Shannon ~ Region/Former_landuse, data = BacIndAgriGenera.rarefied2.merged_replica.df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33604 -0.11658  0.02269  0.09271  0.40594 
## 
## Coefficients:
##                                         Estimate Std. Error t value Pr(>|t|)
## (Intercept)                              4.44821    0.03803 116.951  < 2e-16
## RegionScotland                          -0.18271    0.05379  -3.397  0.00126
## RegionEngland:Former_landuseIndustrial   0.06363    0.05379   1.183  0.24181
## RegionScotland:Former_landuseIndustrial  0.29169    0.05379   5.423 1.29e-06
##                                            
## (Intercept)                             ***
## RegionScotland                          ** 
## RegionEngland:Former_landuseIndustrial     
## RegionScotland:Former_landuseIndustrial ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1473 on 56 degrees of freedom
## Multiple R-squared:  0.3782, Adjusted R-squared:  0.3449 
## F-statistic: 11.36 on 3 and 56 DF,  p-value: 6.367e-06

6 Differential Abundance

On merged sites replicas

library(ANCOMBC)

Differential Abundance Tests: pH, conductivity, Region, age, FLU

pH + conductivity + age –> Trend test Region –> Global test / pairwise FLU –> pairwise

6.1 FLU

# Estimate differential abundances
# 1. FLU (Pairwise) bonferroni correction is often applied to control for the increased risk of Type I errors
ancom.FLU = ancombc2(data = BacIndAgriGenera.rarefied2.merged_replica, tax_level = "Genus",
                     fix_formula = "Former_landuse", 
                     group = "Former_landuse", 
                     p_adj_method = "bonferroni", alpha = 0.05, pairwise = TRUE,
                     prv_cut = 0.10, lib_cut = 0, 
                     n_cl = 4, verbose = TRUE)
# alpha: significance cut-off
# prv_cut: filter taxa not present in prv_cut * samples_number
# lib_cut: filter samples with library size threshold
# Extract the results table
ancom.FLU.df <- ancom.FLU$res
colnames(ancom.FLU.df)
##  [1] "taxon"                              "lfc_(Intercept)"                   
##  [3] "lfc_Former_landuseIndustrial"       "se_(Intercept)"                    
##  [5] "se_Former_landuseIndustrial"        "W_(Intercept)"                     
##  [7] "W_Former_landuseIndustrial"         "p_(Intercept)"                     
##  [9] "p_Former_landuseIndustrial"         "q_(Intercept)"                     
## [11] "q_Former_landuseIndustrial"         "diff_(Intercept)"                  
## [13] "diff_Former_landuseIndustrial"      "passed_ss_(Intercept)"             
## [15] "passed_ss_Former_landuseIndustrial"
# Select only columns that we need 
ancom.FLU.df <- ancom.FLU.df %>% 
  select(taxon, contains("Former_landuse"))
# Select differentially abundant taxa
dif.FLU.df <- ancom.FLU.df %>% 
  filter(diff_Former_landuseIndustrial & passed_ss_Former_landuseIndustrial) %>% 
  arrange(desc(lfc_Former_landuseIndustrial))
# Filter out values that are not 2 fold change in both positive and negative direction
dif.FLU.df <- dif.FLU.df %>%
  filter(lfc_Former_landuseIndustrial >= 1 | lfc_Former_landuseIndustrial <= -1)

# Create the 'direct' column to categorize LFC
dif.FLU.df <- dif.FLU.df %>%
  mutate(
    direct = ifelse(lfc_Former_landuseIndustrial >= 1, "Positive LFC", "Negative LFC")
  )

# Ensure direct factorized 
dif.FLU.df$direct = 
  factor(dif.FLU.df$direct, levels = c("Positive LFC", "Negative LFC"))
dif.FLU.df$taxon = 
  factor(dif.FLU.df$taxon, levels = dif.FLU.df$taxon)
dif.FLU.df
##                    taxon lfc_Former_landuseIndustrial
## 1         Granulicella_B                     1.852967
## 2                 VAYN01                     1.294547
## 3  Sphingomonas_L_486704                     1.136363
## 4           Georgfuchsia                    -1.037880
## 5                Bog-159                    -1.056645
## 6              Gp13-AA74                    -1.142105
## 7         Clostridium_AD                    -1.150575
## 8                   AV80                    -1.333078
## 9             Bacillus_A                    -1.345840
## 10          Sporosarcina                    -1.389808
## 11           Bacillus_BK                    -1.465648
##    se_Former_landuseIndustrial W_Former_landuseIndustrial
## 1                    0.1366279                  13.562139
## 2                    0.2358824                   5.488104
## 3                    0.1619448                   7.016978
## 4                    0.1380735                  -7.516871
## 5                    0.2427298                  -4.353176
## 6                    0.1415131                  -8.070667
## 7                    0.2120892                  -5.424958
## 8                    0.1850792                  -7.202746
## 9                    0.1663562                  -8.090111
## 10                   0.2340446                  -5.938220
## 11                   0.1507787                  -9.720521
##    p_Former_landuseIndustrial q_Former_landuseIndustrial
## 1                9.970687e-06               4.127864e-03
## 2                1.063179e-06               4.401563e-04
## 3                6.082429e-06               2.518126e-03
## 4                5.881795e-07               2.435063e-04
## 5                5.884443e-05               2.436160e-02
## 6                3.237361e-07               1.340267e-04
## 7                1.230925e-06               5.096030e-04
## 8                2.127503e-09               8.807864e-07
## 9                1.914571e-08               7.926325e-06
## 10               1.157035e-06               4.790125e-04
## 11               2.578196e-05               1.067373e-02
##    diff_Former_landuseIndustrial passed_ss_Former_landuseIndustrial
## 1                           TRUE                               TRUE
## 2                           TRUE                               TRUE
## 3                           TRUE                               TRUE
## 4                           TRUE                               TRUE
## 5                           TRUE                               TRUE
## 6                           TRUE                               TRUE
## 7                           TRUE                               TRUE
## 8                           TRUE                               TRUE
## 9                           TRUE                               TRUE
## 10                          TRUE                               TRUE
## 11                          TRUE                               TRUE
##          direct
## 1  Positive LFC
## 2  Positive LFC
## 3  Positive LFC
## 4  Negative LFC
## 5  Negative LFC
## 6  Negative LFC
## 7  Negative LFC
## 8  Negative LFC
## 9  Negative LFC
## 10 Negative LFC
## 11 Negative LFC
# Make bar plot
dif.FLU.plot <- dif.FLU.df %>%
    ggplot(aes(x = taxon, y = lfc_Former_landuseIndustrial, fill = direct)) + 
    geom_bar(stat = "identity", width = 0.7, color = "black", 
             position = position_dodge(width = 0.4)) +
    geom_errorbar(aes(ymin = lfc_Former_landuseIndustrial - se_Former_landuseIndustrial, 
                      ymax = lfc_Former_landuseIndustrial + se_Former_landuseIndustrial), 
                  width = 0.2, position = position_dodge(0.05), 
                  color = "black") + 
    labs(x = "Genus", y = "Log fold change") +  
    ggtitle(label = "Differentially abundant taxa", 
            subtitle="Prior land use: Industrial vs Agricultural") +
    scale_fill_discrete(name = NULL) +
    scale_color_discrete(name = NULL) +
    theme_bw() + 
    theme(panel.grid.minor.y = element_blank(),
          axis.text.x = element_text(size = 6, angle = 60, hjust = 1))

# Convert to plotly 
dif.FLU.plot_plotly <- ggplotly(dif.FLU.plot, tooltip = "text")
dif.FLU.plot_plotly
# Save the plot as an HTML file
htmlwidgets::saveWidget(dif.FLU.plot_plotly, "results/plotly/differential_abundant_taxa_FLU.html")
# box plot function
diff_abund.box_plot <- function(taxon, env_factor, title) {
  taxon_data <- psmelt(BacIndAgriGenera.rarefied) %>%
    filter(Genus == taxon)
  
  plot <- ggplot(taxon_data, aes_string(x = env_factor, y = "Abundance", fill = env_factor)) +
    geom_violin() +
    geom_jitter(width = 0.2, alpha = 0.5) +
    labs(title = title, x = env_factor, y = "Abundance") +
    theme_minimal() +
    theme(legend.position = "none")
  
  return(plot)
}
# Extract the first positive and negative taxa
first_positive_taxa.FLU <- dif.FLU.df %>%
  filter(direct == "Positive LFC") %>%
  slice(1) %>%
  pull(taxon)

first_negative_taxa.FLU <- dif.FLU.df %>%
  filter(direct == "Negative LFC") %>%
  slice(1) %>%
  pull(taxon)


# Generate and plot the box plots
positive_plot.FLU <- diff_abund.box_plot(first_positive_taxa.FLU, "Former_landuse",
                                   paste("Abundance of", first_positive_taxa.FLU, "by Former Land Use"))
negative_plot.FLU <- diff_abund.box_plot(first_negative_taxa.FLU, "Former_landuse",
                                   paste("Abundance of", first_negative_taxa.FLU, "by Former Land Use"))


print(positive_plot.FLU)
print(negative_plot.FLU)

6.2 Region

# 2. Region
ancom.region = ancombc2(data = BacIndAgriGenera.rarefied2.merged_replica, tax_level = "Genus",
                     fix_formula = "Region", 
                     group = "Region", 
                     p_adj_method = "bonferroni", alpha = 0.05, pairwise = TRUE,
                     prv_cut = 0.10, lib_cut = 0, 
                     n_cl = 4, verbose = TRUE)
# Extract the results table
ancom.region.df <- ancom.region$res
# Select only columns that we need 
ancom.region.df <- ancom.region.df %>% 
  select(taxon, contains("Region"))
# Select differentially abundant taxa
dif.region.df <- ancom.region.df %>% 
  filter(diff_RegionScotland & passed_ss_RegionScotland) %>% 
  arrange(desc(lfc_RegionScotland))
sum(dif.region.df$lfc_RegionScotland > 1)
## [1] 17
# Filter out values that are not 2 fold change in both positive and negative direction
dif.region.df <- dif.region.df %>%
  filter(lfc_RegionScotland >= 1 | lfc_RegionScotland <= -1)

# Create the 'direct' column to categorize LFC
dif.region.df <- dif.region.df %>%
  mutate(
    direct = ifelse(lfc_RegionScotland >= 1, "Positive LFC", "Negative LFC")
  )
# Ensure direct factorized 
dif.region.df$direct = 
  factor(dif.region.df$direct, levels = c("Positive LFC", "Negative LFC"))
dif.region.df$taxon = 
  factor(dif.region.df$taxon, levels = dif.region.df$taxon)
dif.region.df
##                    taxon lfc_RegionScotland se_RegionScotland W_RegionScotland
## 1                 VBOC01           1.642275         0.1719702         9.549768
## 2              Palsa-189           1.511405         0.2475229         6.106121
## 3            Nemorincola           1.509296         0.1429288        10.559774
## 4  Flaviaesturariibacter           1.409598         0.2144276         6.573773
## 5          Pelomicrobium           1.391422         0.2058517         6.759339
## 6        2013-40CM-41-45           1.381802         0.1425396         9.694161
## 7        SCGC-AG-212-J23           1.264733         0.1737948         7.277162
## 8            Acidoferrum           1.208742         0.2791036         4.330801
## 9              Gp1-AA122           1.195852         0.2403979         4.974470
## 10              CADDZI01           1.182387         0.1965815         6.014746
## 11               UBA4416           1.178866         0.1712315         6.884633
## 12        Hypericibacter           1.158846         0.1634360         7.090519
## 13          Chitinophaga           1.156916         0.1527348         7.574666
## 14             Gp1-AA133           1.152696         0.2020446         5.705155
## 15             Palsa-187           1.114763         0.1918483         5.810648
## 16          RBG-16-71-46           1.047354         0.2228241         4.700364
## 17                UBA923           1.022077         0.1429087         7.151960
## 18          Actinoplanes          -1.045388         0.1957841        -5.339495
## 19                VFJQ01          -1.048863         0.2197708        -4.772531
## 20        Methyloligella          -1.096814         0.1415061        -7.751000
## 21            Luteitalea          -1.129800         0.2329756        -4.849433
## 22                SHVA01          -1.143939         0.2414934        -4.736937
## 23               HRBIN40          -1.168912         0.2587924        -4.516792
## 24 Streptomyces_G_399870          -1.194410         0.1707691        -6.994297
## 25          Povalibacter          -1.249674         0.2262127        -5.524331
## 26      Parviterribacter          -1.351487         0.1549108        -8.724291
## 27 Microlunatus_A_391020          -1.353893         0.2000433        -6.768001
## 28       Solirubrobacter          -1.405081         0.2846144        -4.936789
## 29                 AC-51          -1.465112         0.2722885        -5.380735
## 30                VFJN01          -1.559506         0.2350581        -6.634555
## 31 Nocardioides_A_392796          -1.610780         0.2123895        -7.584086
## 32             Kribbella          -1.665643         0.2004193        -8.310791
## 33              JAAAPF01          -1.705734         0.1759461        -9.694642
## 34            GMQP-bins7          -1.734388         0.2663932        -6.510632
## 35          Ohtaekwangia          -1.861354         0.2649990        -7.024001
##    p_RegionScotland q_RegionScotland diff_RegionScotland
## 1      1.805978e-08     7.476749e-06                TRUE
## 2      1.605755e-07     6.647824e-05                TRUE
## 3      1.678942e-10     6.950821e-08                TRUE
## 4      1.197430e-07     4.957361e-05                TRUE
## 5      8.640151e-09     3.577023e-06                TRUE
## 6      6.912439e-05     2.861750e-02                TRUE
## 7      6.712846e-09     2.779118e-06                TRUE
## 8      6.489750e-05     2.686757e-02                TRUE
## 9      7.015954e-06     2.904605e-03                TRUE
## 10     1.093486e-05     4.527032e-03                TRUE
## 11     2.134447e-07     8.836609e-05                TRUE
## 12     3.681885e-06     1.524300e-03                TRUE
## 13     4.808854e-10     1.990865e-07                TRUE
## 14     6.263632e-07     2.593143e-04                TRUE
## 15     3.424948e-07     1.417928e-04                TRUE
## 16     2.298778e-05     9.516941e-03                TRUE
## 17     4.919953e-06     2.036860e-03                TRUE
## 18     4.921712e-06     2.037589e-03                TRUE
## 19     1.741330e-05     7.209107e-03                TRUE
## 20     1.114786e-04     4.615212e-02                TRUE
## 21     2.690206e-05     1.113745e-02                TRUE
## 22     1.663815e-05     6.888195e-03                TRUE
## 23     3.542722e-05     1.466687e-02                TRUE
## 24     2.996486e-09     1.240545e-06                TRUE
## 25     1.324347e-06     5.482796e-04                TRUE
## 26     1.100480e-05     4.555987e-03                TRUE
## 27     5.075079e-08     2.101083e-05                TRUE
## 28     8.904792e-06     3.686584e-03                TRUE
## 29     2.884837e-06     1.194323e-03                TRUE
## 30     4.362516e-08     1.806081e-05                TRUE
## 31     3.056924e-10     1.265567e-07                TRUE
## 32     4.121611e-11     1.706347e-08                TRUE
## 33     2.623744e-05     1.086230e-02                TRUE
## 34     1.935353e-08     8.012362e-06                TRUE
## 35     4.987584e-09     2.064860e-06                TRUE
##    passed_ss_RegionScotland       direct
## 1                      TRUE Positive LFC
## 2                      TRUE Positive LFC
## 3                      TRUE Positive LFC
## 4                      TRUE Positive LFC
## 5                      TRUE Positive LFC
## 6                      TRUE Positive LFC
## 7                      TRUE Positive LFC
## 8                      TRUE Positive LFC
## 9                      TRUE Positive LFC
## 10                     TRUE Positive LFC
## 11                     TRUE Positive LFC
## 12                     TRUE Positive LFC
## 13                     TRUE Positive LFC
## 14                     TRUE Positive LFC
## 15                     TRUE Positive LFC
## 16                     TRUE Positive LFC
## 17                     TRUE Positive LFC
## 18                     TRUE Negative LFC
## 19                     TRUE Negative LFC
## 20                     TRUE Negative LFC
## 21                     TRUE Negative LFC
## 22                     TRUE Negative LFC
## 23                     TRUE Negative LFC
## 24                     TRUE Negative LFC
## 25                     TRUE Negative LFC
## 26                     TRUE Negative LFC
## 27                     TRUE Negative LFC
## 28                     TRUE Negative LFC
## 29                     TRUE Negative LFC
## 30                     TRUE Negative LFC
## 31                     TRUE Negative LFC
## 32                     TRUE Negative LFC
## 33                     TRUE Negative LFC
## 34                     TRUE Negative LFC
## 35                     TRUE Negative LFC
# Make bar plot
dif.region.plot <- dif.region.df %>%
    ggplot(aes(x = taxon, y = lfc_RegionScotland, fill = direct)) + 
    geom_bar(stat = "identity", width = 0.7, color = "black", 
             position = position_dodge(width = 0.4)) +
    geom_errorbar(aes(ymin = lfc_RegionScotland - se_RegionScotland, 
                      ymax = lfc_RegionScotland + se_RegionScotland), 
                  width = 0.2, position = position_dodge(0.05), 
                  color = "black") + 
    labs(x = "Genus", y = "Log fold change") +  
    ggtitle(label = "Differentially abundant taxa", 
            subtitle="Region: Scotland vs England") +
    scale_fill_discrete(name = NULL) +
    scale_color_discrete(name = NULL) +
    theme_bw() + 
    theme(panel.grid.minor.y = element_blank(),
          axis.text.x = element_text(size = 6, angle = 60, hjust = 1))

# Convert to plotly 
dif.region.plot_plotly <- ggplotly(dif.region.plot, tooltip = "text")
dif.region.plot_plotly
# Save the plot as an HTML file
htmlwidgets::saveWidget(dif.region.plot_plotly, "results/plotly/differential_abundant_taxa_Region.html")
# Extract the first positive and negative taxa
first_positive_taxa.region <- dif.region.df %>%
  filter(direct == "Positive LFC") %>%
  slice(1) %>%
  pull(taxon)

first_negative_taxa.region <- dif.region.df %>%
  filter(direct == "Negative LFC") %>%
  slice(1) %>%
  pull(taxon)


# Generate and plot the box plots
positive_plot.region <- diff_abund.box_plot(first_positive_taxa.region, "Region",
                                   paste("Abundance of", first_positive_taxa.region, "by Region"))
negative_plot.region <- diff_abund.box_plot(first_negative_taxa.region, "Region",
                                   paste("Abundance of", first_negative_taxa.region, "by Region"))


print(positive_plot.region)
print(negative_plot.region)

7 Ordination

7.1 NMDS

# Compute unweighted UniFrac distance
unifrac_dist <- UniFrac(BacIndAgriGenera.rarefied2.merged_replica, weighted = FALSE)

# Perform NMDS ordination
ordination.NMDS <- ordinate(BacIndAgriGenera.rarefied2.merged_replica, method = "NMDS", distance = unifrac_dist)

## For 3d plotting
# NMDS ordination
ordination.NMDS.3d <- metaMDS(unifrac_dist, k=3)

# Extract NMDS samples scores (the position of each sample in the ordination space)
nmds_scores <- scores(ordination.NMDS.3d, display="sites")

# Convert to data frame and add sample metadata
nmds_data <- as.data.frame(nmds_scores)

###Region

# 2d plot
ordination_NMDS.Region <- plot_ordination(BacIndAgriGenera.rarefied2.merged_replica, ordination.NMDS, type="samples", color="Region") + 
  ggtitle("samples ordination:", 
          subtitle="NMDS with unweighted UniFrac distances")+ 
  geom_point(size=3) + stat_ellipse() +
  theme_bw() +
  theme(legend.position = "bottom") +
  guides(color=guide_legend(ncol=3))

print(ordination_NMDS.Region)

nmds_data$Region <- sample_data(BacIndAgriGenera.rarefied2.merged_replica)$Region

# 3D plot
ordination_NMDS.Region.3d <- plot_ly(nmds_data, x = ~NMDS1, y = ~NMDS2, z = ~NMDS3, 
               color = ~Region, colors = c("red", "blue"), 
               type = 'scatter3d', mode = 'markers',
               marker = list(size = 5))

ordination_NMDS.Region.3d <- ordination_NMDS.Region.3d %>% layout(title = "NMDS with Unweighted UniFrac Distances",
                      scene = list(xaxis = list(title = 'NMDS1'),
                                   yaxis = list(title = 'NMDS2'),
                                   zaxis = list(title = 'NMDS3')))

ordination_NMDS.Region.3d

###FLU

ordination_NMDS.FLU <- plot_ordination(BacIndAgriGenera.rarefied2.merged_replica, ordination.NMDS, type="samples", color="Former_landuse") + 
  ggtitle("samples ordination:", 
          subtitle="NMDS with unweighted UniFrac distances")+ 
  geom_point(size=3) + stat_ellipse() +
  theme_bw() +
  theme(legend.position = "bottom") +
  guides(color=guide_legend(ncol=3))

print(ordination_NMDS.FLU)

nmds_data$Former_landuse <- sample_data(BacIndAgriGenera.rarefied2.merged_replica)$Former_landuse

# 3D plot
ordination_NMDS.FLU.3d <- plot_ly(nmds_data, x = ~NMDS1, y = ~NMDS2, z = ~NMDS3, 
               color = ~Former_landuse, colors = c("red", "blue"), 
               type = 'scatter3d', mode = 'markers',
               marker = list(size = 5))

ordination_NMDS.FLU.3d <- ordination_NMDS.FLU.3d %>% layout(title = "NMDS with Unweighted UniFrac Distances",
                      scene = list(xaxis = list(title = 'NMDS1'),
                                   yaxis = list(title = 'NMDS2'),
                                   zaxis = list(title = 'NMDS3')))

ordination_NMDS.FLU.3d

7.2 RDA

# Convert BacIndAgriGenera.rarefied2.merged_replica to data frame for abundance and environmental data
abund_data <- as.data.frame(t(otu_table(BacIndAgriGenera.rarefied2.merged_replica)))
env_data <- data.frame(as(sample_data(BacIndAgriGenera.rarefied2.merged_replica), "data.frame"))
# Select variables 
env_data <- env_data[, c("Region", "Former_landuse", "Woodland_age", "pH", "Conductivity")]

# Make sure pH, age, and conductivity are numerical
env_data$Woodland_age <- as.numeric(as.character(env_data$Woodland_age))
env_data$Conductivity <- as.numeric(as.character(env_data$Conductivity))
env_data$pH <- as.numeric(as.character(env_data$pH))

# Remove NA values
env_data <- na.omit(env_data)
# Factorize categorical columns
env_data$Region <- as.factor(env_data$Region)
env_data$Former_landuse <- as.factor(env_data$Former_landuse)


# Select row names in abund_data that match env_data row names
abund_data <- abund_data[rownames(abund_data) %in% rownames(env_data), ]

# Perform RDA, constraining all env_data variables
ordination.RDA <- rda(abund_data ~ pH + Conductivity, data = env_data)
RDA.R2adj <- RsquareAdj(ordination.RDA)$adj.r.squared
RDA.R2adj
## [1] 0.356577
# # Enable the r-universe repo
# options(repos = c(
#     fawda123 = 'https://fawda123.r-universe.dev',
#     CRAN = 'https://cloud.r-project.org'))
# 
# # Install ggord
# install.packages('ggord')
library(ggord)
ggord(ordination.RDA, env_data$Region, ext=0.8, txt=4, grp_title="Region", xlims=c(-1,1), ylims=c(-1,1), alpha_el=0.6, size=1.5, poly=FALSE,
      vectyp="longdash", obslab=FALSE, ptslab=FALSE)  

ggord(ordination.RDA, env_data$Former_landuse, ext=0.8, txt=4, grp_title="Former_landuse", xlims=c(-1,1), ylims=c(-1,1), poly=FALSE, size=1.5,
      vectyp="longdash", obslab=FALSE, ptslab=FALSE)  

7.3 cCA

# Perform CCA
ordination.cCA <- cca(abund_data ~ pH + Conductivity, data = env_data)
cCA.R2adj <- RsquareAdj(ordination.cCA)$adj.r.squared
cCA.R2adj
## [1] 0.2424402
ggord(ordination.cCA, env_data$Region, ext=0.8, txt=4, grp_title="Region", xlims=c(-1,1), ylims=c(-1,1), alpha_el=0.6, size=1.5,
      vectyp="longdash", obslab=FALSE)  

ggord(ordination.cCA, env_data$Former_landuse, ext=0.8, txt=4, grp_title="Former_landuse", xlims=c(-1,1), ylims=c(-1,1), alpha_el=0.6, size=1.5,
      vectyp="longdash", obslab=FALSE)  

8 Beta Diversity (nested design)

# Re transpose abund_data
otu_table_matrix <- as.matrix(t(abund_data))

# Create a new phyloseq object with the filtered samples (based on env_data, no NaNs)
BacIndAgriGenera.rarefied2.merged.filtered <- phyloseq(
  otu_table(otu_table_matrix, taxa_are_rows = TRUE),
  sample_data(env_data),
  tax_table(BacIndAgriGenera.rarefied2.merged_replica),
  phy_tree(BacIndAgriGenera.rarefied2.merged_replica)
)
# Calculate beta diversity metrics
weighted_unifrac_dist.filtered <- UniFrac(BacIndAgriGenera.rarefied2.merged.filtered, weighted = TRUE)
unweighted_unifrac_dist.filtered <- UniFrac(BacIndAgriGenera.rarefied2.merged.filtered, weighted = FALSE)
bray_dist.filtered <- distance(BacIndAgriGenera.rarefied2.merged.filtered, method="bray")
adjus_p_value <- function(adonis_result){
  # Extract p-values 
  p_values <- adonis_result$`Pr(>F)`

  # Adjust the p-values
  p_adjusted <- p.adjust(p_values)  

  # Add the adjusted p-values to adonis2 table
  adonis_result$`Pr(>F)_adjusted` <- c(p_adjusted)
  return(adonis_result)
}
library(scales)  

nested_model.pie <- function(nested_model, title){
  # Extract the explained variance (R2) 
  explained_variance <- nested_model$R2
  
  # Remove the NA values associated with the rest of table info
  p_values <- nested_model$`Pr(>F)`
  explained_variance <- explained_variance[!is.na(p_values)]
  
  # Names of the factors
  factors <- rownames(nested_model)[1:length(explained_variance)]
  
  # Create a data frame for plotting
  plot_data <- data.frame(
    Factor = factors,
    ExplainedVariance = explained_variance
  )

  # Ensure the levels of Factor are in the same order as the explained variance
  plot_data$Factor <- factor(plot_data$Factor, levels = plot_data$Factor)
  
  # Custom legends
  plot_data$LegendLabel <- paste0(plot_data$Factor, 
                                  " (", percent(plot_data$ExplainedVariance), ")")

  # Create the pie chart
  pie_chart <- ggplot(plot_data, aes(x = "", y = ExplainedVariance, fill = Factor)) +
    geom_bar(width = 1, stat = "identity") +
    coord_polar(theta = "y") +
    scale_y_continuous(labels = percent_format()) +
    theme_void() +
    theme(legend.position = "right") +
    labs(title = title, fill = "Factor") +
    scale_fill_manual(values = scales::brewer_pal(palette = "Set3")(length(plot_data$Factor)),
                      labels = plot_data$LegendLabel) 

  print(pie_chart)

}

pH / conductivity / FLU / woodland age / region

8.1 Weighted UniFrac

nested_model.weighted_unifrac <- adonis2(weighted_unifrac_dist.filtered ~ Region/Former_landuse 
                                         + pH + Woodland_age,
        data = env_data, strata = env_data$Region, permutations=1000)
nested_model.weighted_unifrac
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 1000
## 
## adonis2(formula = weighted_unifrac_dist.filtered ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
##                       Df SumOfSqs      R2       F   Pr(>F)    
## Region                 1  0.85916 0.31705 46.3139 0.000999 ***
## pH                     1  0.58442 0.21567 31.5038 0.000999 ***
## Woodland_age           1  0.03752 0.01385  2.0225 0.088911 .  
## Region:Former_landuse  2  0.24555 0.09061  6.6182 0.000999 ***
## Residual              53  0.98319 0.36282                     
## Total                 58  2.70984 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adjus_p_value(nested_model.weighted_unifrac)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 1000
## 
## adonis2(formula = weighted_unifrac_dist.filtered ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
##                       Df SumOfSqs      R2       F   Pr(>F) Pr(>F)_adjusted   
## Region                 1  0.85916 0.31705 46.3139 0.000999        0.003996 **
## pH                     1  0.58442 0.21567 31.5038 0.000999        0.003996 **
## Woodland_age           1  0.03752 0.01385  2.0225 0.088911        0.088911 . 
## Region:Former_landuse  2  0.24555 0.09061  6.6182 0.000999        0.003996 **
## Residual              53  0.98319 0.36282                                    
## Total                 58  2.70984 1.00000                                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nested_model.pie(nested_model.weighted_unifrac, "Explained Variance of weighted UniFrac data by environmental factors")

8.2 Unweighted UniFrac

nested_model.unweighted_unifrac <- adonis2(unweighted_unifrac_dist.filtered ~ Region/Former_landuse 
                                         + pH + Woodland_age,
        data = env_data, strata = env_data$Region, permutations=1000)
nested_model.unweighted_unifrac
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 1000
## 
## adonis2(formula = unweighted_unifrac_dist.filtered ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
##                       Df SumOfSqs      R2       F   Pr(>F)    
## Region                 1   0.7242 0.19127 17.4608 0.000999 ***
## pH                     1   0.4694 0.12396 11.3164 0.000999 ***
## Woodland_age           1   0.0719 0.01900  1.7345 0.049950 *  
## Region:Former_landuse  2   0.3226 0.08520  3.8889 0.000999 ***
## Residual              53   2.1982 0.58057                     
## Total                 58   3.7863 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adjus_p_value(nested_model.unweighted_unifrac)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 1000
## 
## adonis2(formula = unweighted_unifrac_dist.filtered ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
##                       Df SumOfSqs      R2       F   Pr(>F) Pr(>F)_adjusted   
## Region                 1   0.7242 0.19127 17.4608 0.000999        0.003996 **
## pH                     1   0.4694 0.12396 11.3164 0.000999        0.003996 **
## Woodland_age           1   0.0719 0.01900  1.7345 0.049950        0.049950 * 
## Region:Former_landuse  2   0.3226 0.08520  3.8889 0.000999        0.003996 **
## Residual              53   2.1982 0.58057                                    
## Total                 58   3.7863 1.00000                                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nested_model.pie(nested_model.unweighted_unifrac, "Explained Variance of unweighted UniFrac data by environmental factors")

8.3 Bray curtis

nested_model.bray <- adonis2(bray_dist.filtered ~ Region/Former_landuse 
                                         + pH + Woodland_age,
        data = env_data, strata = env_data$Region, permutations=1000)
nested_model.bray
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 1000
## 
## adonis2(formula = bray_dist.filtered ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
##                       Df SumOfSqs      R2       F   Pr(>F)    
## Region                 1   1.1178 0.22339 28.1794 0.000999 ***
## pH                     1   1.2356 0.24694 31.1505 0.000999 ***
## Woodland_age           1   0.1036 0.02070  2.6115 0.027972 *  
## Region:Former_landuse  2   0.4445 0.08883  5.6028 0.000999 ***
## Residual              53   2.1023 0.42014                     
## Total                 58   5.0038 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adjus_p_value(nested_model.bray)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 1000
## 
## adonis2(formula = bray_dist.filtered ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
##                       Df SumOfSqs      R2       F   Pr(>F) Pr(>F)_adjusted   
## Region                 1   1.1178 0.22339 28.1794 0.000999        0.003996 **
## pH                     1   1.2356 0.24694 31.1505 0.000999        0.003996 **
## Woodland_age           1   0.1036 0.02070  2.6115 0.027972        0.027972 * 
## Region:Former_landuse  2   0.4445 0.08883  5.6028 0.000999        0.003996 **
## Residual              53   2.1023 0.42014                                    
## Total                 58   5.0038 1.00000                                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nested_model.pie(nested_model.bray, "Explained Variance of bray curtis by environmental factors")

8.4 Abundance data

nested_model.abund <- adonis2(abund_data ~ Region/Former_landuse 
                                         + pH + Woodland_age,
        data = env_data, strata = env_data$Region, permutations=1000)
nested_model.abund
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 1000
## 
## adonis2(formula = abund_data ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
##                       Df SumOfSqs      R2       F   Pr(>F)    
## Region                 1   1.1178 0.22339 28.1794 0.000999 ***
## pH                     1   1.2356 0.24694 31.1505 0.000999 ***
## Woodland_age           1   0.1036 0.02070  2.6115 0.032967 *  
## Region:Former_landuse  2   0.4445 0.08883  5.6028 0.000999 ***
## Residual              53   2.1023 0.42014                     
## Total                 58   5.0038 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adjus_p_value(nested_model.abund)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 1000
## 
## adonis2(formula = abund_data ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
##                       Df SumOfSqs      R2       F   Pr(>F) Pr(>F)_adjusted   
## Region                 1   1.1178 0.22339 28.1794 0.000999        0.003996 **
## pH                     1   1.2356 0.24694 31.1505 0.000999        0.003996 **
## Woodland_age           1   0.1036 0.02070  2.6115 0.032967        0.032967 * 
## Region:Former_landuse  2   0.4445 0.08883  5.6028 0.000999        0.003996 **
## Residual              53   2.1023 0.42014                                    
## Total                 58   5.0038 1.00000                                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nested_model.pie(nested_model.abund, "Explained Variance of abundance data by environmental factors")